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A computational fluid model is developed to study waves and instabilities. A new technique 
involving initial perturbations in configuration space have been implemented to excite the plasma 
waves; i.e. the perturbations acting similar to a random velocity distribution in particle in cell (PIC) 
codes. This forms a new powerful tool for investigation of many waves arising in both thermal and 
cold plasmas and as such will allow investigation of problems demanding scales and resolution not yet 
possible by PIC codes. The model predicts Langmuir waves, two stream instabilities, nonlinear wave- 
wave interaction, and the Debye screening effects. The agreement between theory and simulation 
where analytic results are available are excellent. 

I. INTRODUCTION 

Due to the fundamental importance of the waves and instabilities in plasma and hydrodynamics investigations, 
computational researchers have devoted great efforts in developing appropriate tools. One of the main challenges after 
developing numerically stable algorithms in fluid models has been generation of the waves in the linear, nonlinear as 
well as unstable modes; i.e. waves which preserve analytic dispersion relations 1 [2]. Furthermore extending the case 
of hydrodynamics to that of MHD and or plasma physics one deals with waves with considerably more complicated 
propagation characteristics than the hydrodynamics cases treated by those authors; i.e. dispersion, polarization, 
oblique propagations, etc. 

The main problems in generating a wave spectrum from small amplitude disturbances in fluid equations are: (1) the 
highly nonlinear nature of those equations; (2) the lack of an initial thermal velocity distribution. The first problem 
could cause any small amplitude configuration space disturbance to grow to very large amplitudes in relatively short 
times and result in wave breaking and non-propagation. Also when there does exist a thermal distribution, there 
are always a distribution of thermalized particles in phase with most waves; they can therefore excite the allowed 
modes to at least half their thermal level. Therefore in a case without thermal equilibrium, a disturbance of arbitrary 
wavelength cannot strictly speaking apportion its energy to other allowed modes. For example in purely electrostatic 
cases, we know from equilibrium statistical mechanics that when there exist a thermal distribution each mode Ei{k) 
can acquire an energy [3]: 

< E f ' 2> oc kT. (1) 
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1 The waves'propagation characteristics are encoded in the dispersion relations [1] 
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To investigate MHD wave spectra therefore magnetohydrodynamic particle codes have served as powerful tools [4], 
[6], [7], and [8]. For other plasma waves PIC [9] and [10] or hybrid codes [8], [11], and [12] have served as the main 
wave investigation tools; i.e., basically codes which start from thermal equilibrium. In these codes the random particle 
distribution acts like a disturbance in velocity space and configuration space remains unaltered at the beginning of 
each simulation. 

In our case we initiate each simulation by a perturbation in configuration space. Despite the initial shape of the 
perturbation, we observe other allowed modes to develop similar to PIC simulations. We believe that the mesh dis- 
cretization and the finite differencing contribute in the following ways: (i) round of errors alter the initial perturbation 
shape and can drive other wavelength; (ii) as the nonlinear effects grow amplitudes and shorten wavelengths to the 
numerical dissipation and dispersion scale lengths, these effects can act to dampen and initiate the propagation of 
the different modes and prevent indefinite nonlinear growth. These effects can therefore explain the observed wave 
spectra. With this then we can use fluid instead of PIC codes as a convenient alternative to investigate many waves. 

The organization of the paper is as follows: in section II the model is treated analytically; in section III the numerical 
scheme (algorithm, stability and conservation laws) are presented; in section IV the various tests of the model are 
presented (test of the dispersion relation, two stream instability, screening effect and nonlinear harmonic generation). 
At the end a brief summary and conclusion with future direction are presented. 



II. ANALYTICAL TREATMENT 



We focus on the investigation of the high frequency (hf) longitudinal waves; i. e. a frequency domain where ions can 
be safely assumed to form an immobile background (no represents their uniform density). The appropriate equations 
are then Poisson's and the electron fluid equations: 

dn d , . 

m + o-J nv) -^ (2) 

dv d e d 1 dP 
dt dx mdx^ nm dx ' 

—f- = 4 7re (n-n a ). (4) 

Here (p is the self-consistent electric potential, and n, v, P and m represent the electron density, velocity, pressure 
and rest mass respectively. Without any loss of generality this problem is treated in one dimension. These basic 
equations are supplemented by an "equation of state" according to the particular thermodynamic properties of the 
fluid of interest. Here, isothermal equation of state is used: 

P = nT, (5) 

where T is the electron temperature and is assumed to be constant and Boltzmann's constant, k, is assumed to be 
unity. 

The minimum requirement of any computational model lies in its ability to preserve conservation laws; for that 
fluid equations are cast in flux conservative form. Equation (3) in conservative form upon using Eq. (5) in Eq. (3) 
becomes: 
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dv d fl 2 e T 



dt dx\2 



- v 



ip-\ Inn =0. (6) 



Note that the logarithmic term is caused by the electron pressure. Therefore the three equations that form the basis 
of our model are: 

9n d , . 

% + 2 -^ + -lnn)=0, (8) 

dt ox \2 m m J 

d 2 ip 

— =^e{n-n a ). (9) 

We will next derive a dispersion relation for wave propagation using Eqs. (7), (8), and (9). To do this, linearizing Eqs. 
(7), (8), and (9) about a spatially uniform equilibrium (n = ra + Sn, v = Sv and (p — Sip), we obtain the following set: 

dSn d . , „. 

■^-+"0^ = 0, (10) 

^ + ^(-^ + ^n)=0, (11) 
at ox \ m n ) 

^=4neSn. (12) 
Assuming simple plane wave solutions, Eqs. (10), (11), and (12) reduce to the following set of equations: 

— iuoSn + ikn Sv = 0, (13) 

e 1 

-iuj5v + ik( Sip H Sri) = 0, (14) 

m n 

—k 2 Stp = iireSn. (15) 
Eqs. (13), (14), and (15) yield nontrivial solution if the following is obeyed: 

u? = uj 2 p + k 2 v^, (16) 

where 

2 47re 2 n 2 T 

Wp = and = — (17) 



are the electron plasma frequency and the thermal velocity, respectively. 

Studies of Langmuir waves (hf electron waves) are of particular importance. Aside from the applications to real 
experimental situations which will become evident in the application section, they serve as excellent probes for testing 
the validity of the fluid code that we have developed. 
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III. NUMERICAL ALGORITHM 



Our model is simply an intuitive construct based on well-known fluid dynamics and Poisson's equations, geared 
toward plasma physics applications, where many different wave phenomena in dispersive media are of interest. Its 
physical " conceptual basis" can be regarded as a model that treats non-stationary electron wave motion for hf domain 
where u> 3> kvx in linear and nonlinear regions. Besides, it can predict electron wave spectrum more accurately than 
"particle in cell simulation" as here we expect less numerical noise. 



A. Normalization 



In these calculations we use the following normalizations: 



x v n eip , . 

u p t^t, >x, >v, ► n, — -> <p, (18) 

r D v T n T 

where 

r l = Y~~2 — ( 19 ) 

is the electron Debye length. Using these definitions, Eqs. (2), (4), (6), and (16) can now be rewritten as follows: 

m + &<"") = °' (20) 



dt dx \2 

d 2 <p 



8x*~ n - l > (22) 



lo 2 = 1 + k 2 . (23) 



It is already mentioned, logarithmic term in Eq. (21) is caused by the electron pressure.. Thus the code has the 
flexibility of being easily converted to the case when electron pressure is negligible. 



B. The Numerical Scheme 



Next we shall describe the numerical scheme. The steps of the scheme are summarized in Table I. A Lax-Wendroff 
method is used to push n and v, while a poisson solver at the end of each step updates the electric potential. 

The grid spacing and time step are denoted by A and At respectively. The fluid velocity and density are known 
at integer time step I. To complete the initial conditions, ip is computed at the same time step (I) by the help of a 
Poisson solver that is based on tridiagonal matrix method. Then n and v are pushed from I to I + 1/2 as the auxiliary 
step of the Lax-Wendroff scheme using Eqs. (20) and (21) (please refer to item 3 of the Table I). Then again ip is 
computed in the auxiliary step (/ + 1/2) using the value of n in the mentioned step. Having known n, v, and ip at the 
time step 1 + 1/2, we push n and v all the way to time step I + 1 as the main step of the Lax-Wendroff scheme in Eqs. 
(20) and (21) (items 5 and 6 in Table I). The electric potential <p is then computed at the time step I + 1 using n l+1 . 
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TABLE I 

Numerical Algorithm of the Fluid Model for Plasma Waves 



Initially we have: n l m , v l m 

1. Compute electric potential, ip l m , using Poisson solver. 

2. Compute fluxes in continuity and momentum equation in main step: 

(fn) m = n m v mi 

3. Push velocity and density half a time step: 



J+l/2 
m+1/2 

, '+1/2 
y m+l/2 ' 



\{ n m+l + n m) 2A [(/«)' 



m+l mj ? 



-(v 



i 

m+l 



V 1 ) 



At 
2A 



[CM 



m+l 



(/,) 



4. Compute electric potential in half step, <Pm+i/2' usm g n m+i/2' 

5. Compute fluxes in continuity and momentum equations in half step: 



(f V +1/2 - n l+1/2 J +1/2 

If \«+l/2 
\Jv) m+ i/2 



m+1/2 m+1/2' 
1( 1+1/2 yj 1+1/2 +lnn '+l/2 
2^+1/2^ ^m+1/2 + iIln m+l/2- 



6. Push the velocity and density another half a time step: 



J+i 



At 
A 

At 
A 



Un^m+1/2 



(/r 



(/«) 



1 + 1/2 
m+1/2 



J + l/2 
'm-1/2 

(f V + l/2 ' 
\Jv) m _i/2 



7. Compute electric potential in the main step, ip 1 ^ 1 , using n 1 ^ 1 . 



C. Conservation Laws 

Equations (19), (20) are in conservative form, and we demand that the corresponding difference equations to be 
equally conservative. More specifically, we expect finite difference scheme to conserve the mass density (J^^ndx), 
momentum and the energy of the system, irrespective of the errors incurred by the finite difference lattice. 

To investigate the conservation laws, in what follows, a method compatible with both the auxiliary and the main 
steps will be presented [14]. That is, Eqs. (20) and (21) are integrated over each space-time cell (m) of area AtA m 
(At = t l+1 -t l ) as follows: 

A+l A+l 

f f dn f f d 

/ dt / dx— = dt dx—(nv), (24) 

U L m dt J tl A dx y 



t l + 1 Q t l + 1 

\ dt dx^- = - / dt I 
Jt' JA m dt J t i J A 



dl i dx-^- ( -v 2 — ip + Inn 
A ox \2 



(25) 

If JA m Ul Jt> 

Here J A denotes integral over the cell labelled by m. Carrying out trivial integration over dt and dx on the left and 
right sides respectively Eqs. (24) and (25) become: 



/ n l+1 dx- / n l dx=- / dt V\m;) TO , 
JA m JA m Jt> 



(26) 



t i+i 

/ v l+1 dx- v l dx = - dt V ( -v 2 - <p + lni 

J A m JA m Jt l „ V 2 



(27) 
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where a stands for the boundaries of every cell (the right and the left). Using 

(28) 



the following equations are thus obtained: 

.1+1 

f #1 (29) 



n l+1 = n l 



f ! + l 

„i+l 



<r=<-\ dt-y \-v 2 -y + \nn\ . (30) 



Summing over cells (m) in the system results in: 

M M 



E Kt 1 -^4) = - E /, (31) 

tx=1 m=l *' a 



M M t l+i / s 

E(^ +1 -^) = E/ d4Eu- 2 -^+ ln -) • ( 32 ) 

m=l m =l"' ti a V /m 

Since finite differences were used in computing all the derivatives, then if one sums over all the grid cells in the system, 
each such quantities will appear twice with opposite signs corresponding to the cell boundaries that are being shared 
between the neighboring cells, and they will thus add up to zero. There can, however, be contributions from the 
walls of the computation box. For the periodic boundary condition the walls contributions gives zero; for other cases 
appropriate boundary conditions are implemented to insure good conservation using guard cells. 



D. Numerical Stability Analysis 

In order to obtain the Courant-Fredricks-Lewy (CFL) condition for the model, the difference equations (obtained 
from the differential equations for the problem by discretizing them) must be considered. We follow the method 
of Potter [14]; i.e. obtain the integration time pusher operator from the difference equations assuming a spatially 
uniform system and solve them in Fourier space and obtain a non-local result. We shall do the stability analysis with 
the pressure term. 

Recall that the differential equations (20), (21) and (22) formed the basis of the model. These equations upon 
linearization, give: 

dSn d . „ ,„„. 
■W + 8-J V °' (33) 

dSv d 

_ + _(-*, + ft,) =0, (34) 

(35, 
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Next using Eqs. (33), (34), and (35), after combining the auxiliary and the main steps of the Lax-Wendroff scheme 
and assuming n, v, and ip to have the form (I refers to the time step and m inside the parenthesis to the grid location 
along x) 



(n\ v l , p l ) = (n z , v\ p l )e l{kmA \ 
we obtain the following integration matrix (cr = fcA/2): 



(36) 





1+1 ( 


V 






V 



1_ 2A|lsin 2 a 



^A* sin cr cos cr- ^cotcr 1 - 4f - sin 2 <r 



At z 



iAAt 



COt (7 





At 



At' / 

2 / 



\ 






V 




w 



(37) 



Thus, according to Von Neumann stability condition the following inequality should be held: 2 

M < i, 



(38) 



where are the eigenvalues of the integration matrix and subscript refer to different eigenvalues (here /i = 1,2,3). 
The value of g^ is then determined by setting the following determinant equal to zero; i.e., 

-2iAt t 

bin v tus u ^ 2 Bm u 

=0 (39) 



1-^Af^a-g 



-2iAt • iAAt , 1 At 2At • 2 

-£±z±i. sin a cos a - cot cr 1 — =± sin a - g 



At 2 
2 



A 2 
4 sin 2 <r 



iAAt 



cotcr 



2 y 



The corresponding solutions for g are simply: 



.91 = 0, 



02,3=1 - -At 2 - 



1 ■ 2 2A ^ -2 , U o of, 4. 2 , 

^j- suT a ± zJ At 2 cos 2 cr ll + sm CT I • 



5i fulfills the inequality (38). For the two other eigenvalues, we have: 



1.92 1 = |53| = 



I A/- f 1 + A S in 2 cr) + At 4 Q + cos 4 <t J ( 1 + sin 2 a 



1/2 



Equation (38) is then obeyed if the following inequality is held: 



At 2 



1 



COS cr 1 



■ sin 2 cr < 1. 



4 ' ) \~ ' A 2 

Since At and A are small values (0 < At <C 1 and < A < 1) the inequality (42) will be satisfied if: 



At 



< 



(40) 



(41) 



(42) 



(43) 



A " VITA 2 ' 

Inequality (43) is exact up to the scheme accuracy, however, taking into account the smallness of At and A the 
following stability condition results: 



2 Q l+1 = dQ l where g = e ! " A *; Von Neumann stability condition holds for oj real. 
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IV. TESTING THE CODE 



As mentioned, we have constructed the one-dimensional version of the code and have tested it by looking at small 
and large amplitude (nonlinear) effects in an initially uniform plasma. In what follows, a review of the results will be 
given. 

A. Dispersion relation 

The most basic requirement of a computational model aside from conservation laws is its ability to predict the 
linear theory; e.g. the waves dispersion relation. The degree to which the analytic dispersion relation is obeyed acts 
as a gauge of the computational model and serves to determine its limitations. 

From Eq. 40, the dispersion relation of the corresponding difference equation is: 

e-"" At sin(wflAt) = y / (At) 2 cos 2 CT(l + ^sin 2 a). 

where oj = u)r f o>j. Comparison of this with the analytic dispersion relation shows that by changing k — ► 
fcsin(fcA)/(fcA) in the analytic case one roughly recovers the above result for At lor <C 1 , fcA <C 1 . The fact 
that uji does not have any k dependence, implies no part of the k space to be more susceptible to numerical instability 
than others 3 . The difference dispersion relation above also indicates that for sin(fcA)/(fcA) — ► 1 the numerical 
dispersion to disappear; i.e. for modes with wavelengths long compared with the grid spacing it should be negligible. 

For the initial perturbations, small fluctuations in the density from a uniform background were implemented. Table 
2. shows three different initial perturbations used in the simulations; i.e. : 



n(x) 


= I + 0.0Isin(fco 


x) 




n(x) 


= 1 + 0.01(- 


-x + 


X 3 ) 








— I 


fa: 


-1 < x < 


n(x) 


= I + 0.01 < 


I 


- X 


< x < 1 











Else where 



The reason for these choices is that the first perturbation maintain harmonics with wave numbers very close to ko 
while the latter two maintain harmonics more uniformly distributed in the k space. The most important reason for 
such choices was to determine the impact of the initial perturbations on the final wave spectra; strictly speaking the 
latter two are expected to give rise to more uniform spectra. The initial velocity profiles corresponding to these three 
profiles are drawn in Figs. 1(a), (b) and (c). These velocity profiles indicate broader and more uniform distribution of 
bulk flow velocities in the latter two; i.e. the volume of phase space available to wave propagations are considerably 
larger. 



3 Many PIC algorithms show uji oc k 2 ; i.e. intense short wavelength noise or instability. 
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FIG. 1. Velocity profile for a) n{x) = 1 + 0.01(— x + x A )e x , b)n(x) = 1 + 0.01 sin(fcoa;) and c) Saw-tooth function 



Given these two facts though, the plots of the power spectra 4 of the modes versus w (their frequency) indicate very 
close agreement in all the cases; i.e. regardless of the initially excited modes and phase velocities, most the allowed 
fc-space tends to get excited. This supports our earlier claim that the discretization procedure and the numerical 
dispersion and dissipation have in effect broadened and stabilized the initial spectrum. 

Finally the plots of the dispersion relation for a system size of 1024A with A = 0.01 are shown in Fig 2 and 3. The 
close agreement between the analytic theory (solid lines) and the model (circles) for wave numbers k as large as 6 
indicate resolution of the modes with wave lengths of the order of grid spacing with negligible numerical dispersion. 
Comparison of these with the corresponding PIC simulations for a system 256A length (Fig. 4) clearly indicate 
resolution of much shorter wavelengths here and considerably less numerical dispersion. This is understandable since 
in the PIC models the finite particle size effects introduce additional numerical dispersion which cause smaller allowed 
Jfc's. 



4 The power spectrum is determined in two steps: First, the spatial FFT is used in a quantity (e.g. E(x,t)) and stored E(fci,t), 
next for each fa temporal FFT is performed on E(fa,t) 
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FIG. 2. Dispersion relation for Langmuir wave. 




FIG. 3. Dispersion relation for Langmuir wave with Doppler effect 




FIG. 4. Dispersion relation for Langmuir wave for a typical PIC simulation 



10 



One last remark about the cases corresponding to Figs. 2 and 3 is that the latter involves the case in which the bulk 
plasma had an initial flow velocity. Fig. 3 not only shows that the doppler shifted waves also obey their respected 
dispersion relation, it also shows how any "resulting" plasma flow could impact those waves. That is if any nonzero 
average flow should arise from the initial perturbations (i.e. if the scheme does not preserve momentum conservation 
) the dispersion relation would be impacted as in Fig. 3. A glance at Fig. 2 though points that there could not 
have been any doppler shift and therefore no net plasma flow must have resulted from the initial perturbations. 
Calculations also showed that (vf) — initially remained so to round off errors throughout the simulation. So these 
plots also probe the momentum to be conserved in the model. 

B. Wave Launching on the Boundary 

In the next example a wave is launched from the boundary and its behavior is followed. Theoretically, recall that 
in an unmagnetizcd plasma and in the linear regime the plasma shields any incoming AC density perturbation whose 
frequency is less than plasma frequency (a> p ). This effect is shown in Fig. 5(b) and Fig. 6. In this example the 
frequency of the applied density perturbation is half of the plasma frequency. The wave is launched at x = — 25Xd- 
The amplitude of the density perturbation has the following range: nonlinear (0.2,1.8) Fig. 5(a) and linear (0.99,1.01) 
Fig. 5(b) 5 . The penetration depth is from x = (—25, —20) in the linear and x = (—25, —15) in the nonlinear case. 
Furthermore as Fig. 6(a) indicates, upon penetration, after one wave period following the first crest (x = —18), the 
second crest steepens with its wavelength decreasing to grid cell scale. 6 In the linear regime though [Fig. 6(b)] no 
steepening can be seen. 



5 In these particular shots the wave trough fall at launch points. 

6 The oscillations are numerical in nature. The model should be modified to include FCT filter [15] to eliminate these spurious 
oscillations. 
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FIG. 5. Non-linear and linear penetration of electric field (both plots are sketched at t=10). a) Nonlinear case b) linear case 

In the other case, with the same initial condition (respect to linear case), we launched a wave whose frequency 
was larger than the plasma frequency( u > u) p ). This time the density perturbation propagated into the plasma with 
its wavelength and amplitude unchanged as it penetrated the plasma. Its behavior also conformed with the analytic 
dispersion relation. The results are shown in Fig. 7. 
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FIG. 6. Density versus the position when the external frequency is half of the plasma frequency. To give a time evolution 
feeling, they are plotted for five different normalized time. 
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1.015 




FIG. 7. Density versus the position when the external frequency is two times of the plasma frequency. To give a time 
evolution feeling, they are plotted for five different normalized time. 



C. Two Stream Instability 

As a more severe test of the code, we treated the two stream instability. Although the instability arises under 
a wide range of beam conditions, we shall consider only the simple case of two countrastreaming uniform beams of 
electrons with the same number density no- The first beam travels in the x direction with drift velocity v d and the 
second beam in the opposite direction with same drift velocity, i.e. the countrastreaming beams have the same speed. 
The dispersion relation is as follows: 

£ 1 £ = 1 (44) 

{kv d - uj f ^ (kv d + uj f K ' 

where lo 2 — 4ire 2 no/m is the same plasma frequency for both beams. One can then obtain the following expression 
for uj 2 : 

aj 2 =uj 2 p + k 2 v 2 ±u Jp (u J 2 p + Ak 2 v 2 ) 1/2 . (45) 

This relationship between lo 2 and k 2 is shown graphically in Fig. 8. It is clear that, there exists a critical wave number 
k c which separates the stable and unstable modes. In fact , for k 2 < k 2 two values of u) are complex, one of which 
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represents a growing wave; i.e. an instability. Moreover, there exists a wave number k m that corresponds to the most 
unstable mode. 




FIG. 8. Representation of relationship between u> 2 and k 2 . 

These effects are examined by the fluid code. In this case the code was generalized to a two countrastreaming fluid 
model. As the two countrastreaming beams emerging from the opposite ends meet half way into the simulation box, 
a growing wavclike disturbance develops. Figs. 9 and 10 show the evolution of this disturbance for the cases with and 
without the pressure terms respectively. In both cases the disturbance grows locally while in the latter it also begins 
to propagate in both directions; i.e. a result of the dispersion due to the pressure term. 




t-8. 


4 
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A 


i 


( 





FIG. 9. Electric field versus the position in absence of pressure. Time is normalized by lu p . 
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t=0.25 






FIG. 10. Electric field versus the position in presence of pressure. Time is normalized by lu v . 

Furthermore, the instability of each mode was investigated using the mode energy discussed in the previous section: 
i.e. 

P(k,t) = \E(k,t)\ 2 (46) 

The time derivative of this function with respect to k is shown in Fig. 11. As expected, there exists a critical wave 
number bellow which unstable modes can grow. Furthermore we observed the the most unstable mode corresponding 
to k = k m as the maximum in the Fig. 11. Also the dynamic evolution of the beam-beam interaction was observed 
as a movie and both the disturbance growth and upstream propagations (when pressure term was included) were 
observed. 
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k 

FIG. 11. dp(k,t)/dt versus k. Cutoff and maximum wave numbers (k c ,k m ) are comparable with theory. 

V. CONCLUSION 

The result of this paper demonstrates that fluid model can be used to investigate any waves predicted by their 
basic set of equation. This can include waves of kinetic nature with and without dispersion with resolution far greater 
than the corresponding PIC codes. It was demonstrated that appropriate initial perturbations coupled with difference 
algorithms of sufficient but not excessive numerical dispersion and dissipation can give rise to wave spectra spanning 
all the allowed k-space. Many areas of plasma and or space research can greatly benefit from these techniques. 
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